#load required packages
library(EasyABC)
library(abc)
library(bayesplot)
library(ggplot2)
library(abcrf)

#load required data
load("sum_stats.RData")

#set up function for ABC, using a modified version of the transmission function in HERAChp.KandlerCrema that uses less memory
transmission_model <- function(parameter){
  source("hpcc/modified_herachp.R")
  set.seed(parameter[1])
  
  samples <- transmission(N = 729, timesteps = 232, mu = 0.03672527, warmUp = 200, top = 142, bias = parameter[2])
  return(c(samples$x, mean(samples$obs.div)))
  
  rm(samples)
}

#set priors
priors <- list(c("unif", -0.2, 0.2))

#run simulations
abc_transmission <- ABC_rejection(model = transmission_model, prior = priors, nb_simul = 100000, use_seed = TRUE, n_cluster = 11)

#run rejection algorithm
sum_stats_obs <- c(sum_stats[1], sum_stats[2])
abc_rej <- abc(target = sum_stats_obs, param = abc_transmission$param, sumstat = abc_transmission$stats, tol = 0.01, method = "rejection")
